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Abst ract 

A generalized formulation for constructing 
second- and higher-order accurate TVD (total vari- 
ation diminishing) schemes Is presented. A given 
scheme Is made TVD by limiting antldlf f usl ve flux 
differences with some nonlinear functions, 
so-called limiters. The general Idea of the formu- 
lation and Its mathematical proof of Harten's TVD 
conditions Is shown by applying the Lax-Wendroff 
method to a scalar nonlinear equation and constant- 
coefficient system of conservation laws. For the 
system of equations, we derive several definitions 
for the argument used In the limiter function and 
present their performance In numerical experiments. 
Then the formulation Is formally extended to the 
nonlinear system of equations. It Is demonstrated 
that use of the present procedure allows easy 
conversion of existing central or upwind, and 
second-or higher-order differencing schemes so as 
to preserve monotonicity and to yield physically 
admlssable solutions. The formulation Is simple 
mathematically as well as numerically; neither 
matrix- vector multiplication nor Rlemann solver 
Is required. Roughly twice as much computational 
effort Is needed as compared to conventional 
scheme. Although the notion of TVD Is based on 
the Initial value problem, application to the 
steady Euler equations of the formulation Is also 
made. Numerical examples Including various ranges 
of problems show both time- and spatial-accuracy 
In comparison with exact solutions. 

Introduction 

Recently several second-order, nonoscl 1 latory 
schemes have been proposed for solving hyperbolic 
system of conservation laws, as a means to obtain 
accurate weak solutions, (hereafter the term 
"second- order accuracy" Is applied only In the 
region excluding point of extrema In fluxes) See, 
for example, van Leer 1 , Harten 2 , Roe 3 , and Osher 
and Chakravarthy^. Sweby 5 presented a unified 
formulation In which these Independently proposed 
schemes boiled down to using different forms of 
flux limiters due originally to van Leer 1 . 

Harten 2 Introduced the notion of TVD (total vari- 
ation diminishing) and gave sufficient conditions 
for a class of 5-polnt, explicit schemes to be 
second-order accurate as well as TVD. It Is easy 
to show that all second- and higher-order schemes 
can always be recast as a combination of some 
first-order, stable scheme with higher-order terms. 
The first-order scheme Is too dissipative, espe- 
cially crude In the regions of high gradients. 
Consequently these higher-order terms are added to 
take away excessive dlffuslvlty, hence appropri- 
ately called antldlf f usl ve terms by Boris and 
Book®, while maintaining stability. It Is well- 
known that conventional second- (or higher) order, 
shock-capturing schemes produce spurious oscilla- 
tions near dl scontlnutles and generally require 
considerable numerical dissipation. Thus, the 
basic schemes must be altered (e.g., by limiting 


flux differences In our case) to meet such require- 
ments as (1) yielding monotonic and sharp repre- 
sentation of jumps, and (2) satisfying the entropy 
condition so as to automatically select physically 
relevant solutions. Van Leer 1 first successfully 
attacked this problem and arrived at his second- 
order scheme. The crux of the matter Is to form a 
nonlinear combination of the underlying first-order 
scheme and the remaining antldlf fusive terms. 

Hence we view that In second- and higher-order TVD 
schemes the antldlf fusive terms basically play 
seemingly uncompromising roles of both Increasing 
accuracy and diminishing the total variation. We 
also remark that all TVD schemes cited above use 
some form of an upwlndlng scheme as the underlying 
first-order scheme. In this paper we apply this 
principle to formulate our TVD scheme with empha- 
sis on nonlinear systems of conservation laws, 
namely Euler equations. This approach appears 
natural and straightforward for extension to 
multispace dimensions. Although rigorous mathe- 
matical procedures are followed to ensure Harten's 
TVD conditions, the actual Implementation for 
solving the multidimensional Euler equations Is 
quite simple and general. In numerical tests we 
even carried over the same difference schemes, 
proved to be TVD In the hyperbolic system, 
directly to the steady Euler equations whose type 
Is not known. The same characteristics as that 
obtained by solving unsteady equations were found. 

Section 2 shows the present formulation for 
the nonlinear scalar equation. The procedure fol- 
lows closely to that of Sweby 5 , especially In 
the TVD proof, but they differ In the definition 
of a basic first-order scheme. The entropy con- 
dition Is also examined. Section 3 shows the 
extension to the one-dimensional system of con- 
servation laws and the corresponding TVD property 
for a linearized system. Finally, a description 
of the generlzatlon to the Euler equations of gas 
dynamics Is given In Section 4. Application to 
steady equations, use of both explicit and 
Implicit schemes, central and upwlndlng differ- 
ences are also Included. 

Seco nd-or der TVD Scheme for Nonlinear 
Scalar Equat ion 

To faclllate understanding of the present 
formulation, first we consider the numerical solu- 
tion of the scalar nonlinear conservation law In 
one space 

u. + f (u) = 0, t > 0, -«> < x < ® 

t x v 

( 2 . 1 ) 

u(x,0) = U Q ( x ) 

For Eq. (2.1) to be of hyperbolic type, the flux 
function f(u) has a real derivative a = df/du 
with characteristic curve described by dx/dt = a. 

Suppose a second- or higher-order (centered 
or one-sided biased) differencing scheme Is given 
(e.g., Eq. (2.5b)), the present formulation makes 
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It TVD- satisfying by the following steps: (1) 

Form a first-order, upwind differencing, 

(Eq. (2.5a)); (2) Carry out subtraction of the 
dlf ferenclngs and obtain the remaining terms, 
called antldlf fuslve terms, (Eq. (2.6)); (3) Limit 
the antldlf fuslve terms, (Eq. (2.7)), using the 
so-called flux-limiter function 3 ;* and (4) Deter- 
mine the bounds of the limiter from Harten's con- 
ditions 2 , as well as specify Its function forms. 


It Is noted that allowing different functions <p + 
and <p- to be associated with f + and f- Is 
essential for the proof of TVD. 

Now, the Integration scheme for Eq. (2.1) can 
be applied. To Illustrate some preliminaries of 
TVD and Harten's conditions, we consider the 
first-order Euler explicit scheme with the use of 
upwind differencing, Eq. (2.5a). Then we have 


The upwlnding Is conveniently achieved by 
splitting a(u) (a Is real function) Into "+" and 
parts. 


Uj + 1 = Uj - X ( Af j 4- A'fj)", \ = At/AX 

(2.8a) 




where 


a = a + a 


a > 0, a“ < 0 


(2.2a) 


(2.2b) 


or 


, + r - - + 

/A A f . 

u" +1 = U " X[- A aV ♦ — i A-U. " 

3 3 \A tlj 3 A Uj 


(2.8b) 


Thus the flux f Is split accordingly by requir- 
ing: 


(1) f . f ♦ f 


(2.3a) 


(11) a + = (f\. a' = (f') u (2.3b) 


This is rewritten In a general form 


U j +1 " U J ’ C J-1 A U J + C j A+U J 


(2.9) 


We now define the total variation of the mesh 
solution u to be 


We note that how the splitting (not unique) 

Is done Is Irrelevant In the present analysis. It 
is needed only if an actual calculation is made. 

We also remark that Roe's flux difference splitt- 
ing 3 can be applied equally well, as will become 
clear later. 

Let us use the notation 

A + ( )j = ( )jrl - ( )j 

(2.4) 

M >J - < >j - < >j-l 

We now approximate f x : 

f = (A + f~ + A~f + ) /Ax , first-order upwlnding 

(2.5a) 

= (A + f + A f ) /2Ax, central- dlf ferenclng 

(2.5b) 

Subtracting Eqs. (2.5a) from (2.5b), we have the 
anti dif fusive term 

A~ ( A*f f - A + f-)/2Ax (2.6) 

The flux differences A+f 1 are limited by assign- 
ing the limiter functions <p+ and <p- respectively, 
and we write 

f = (A + f" + A~f + )/AX t A ~ (<p + A + f + - <p~aV)/2AX 

(2.7) 


*S1nce It is the flux differences, rather 
than fluxes themselves, that are limited. It might 
be more appropriate to call this the flux- 
difference limiter. 


TV( u) = ZIaVI (2.10) 

i 1 31 

A numerical scheme Is said to be TVD If 

TV(u n+1 ) < TV(u n ) (2.11) 

Harten 2 gave sufficient conditions for Eq. (2.9) 
to be TVD. 

Lemma (Harten). Let the coefficients C + and 
C~ In Eq. (2.9) satisfy the Inequalities 

C* > 0 (2.12a) 

C* > 0 (2.12b) 

Cj + c" < 1 (2.12c) 

then the scheme Eq. (2.9) Is TVD. 

For smooth functions f 1 , we have 

— - 1 » a. + 0(A u.) 

A U j 

^ - ♦ 

— J = a. + 0(A u.) 

A u j 

where A + Uj = 0(Ax) 

Hence Eqs. (2.12a) and (2.12b) are satisfied by 
Eq. (2.2b) and the third Inequality, Eq. (2.12c) 
requires 


2 




or 


X a 


= X(aj - aj) < 1 


(2.13a) 


(2.13b) 


which Is Identical to the CFL stability condition. 

Roe 3 proposed another type of splitting In 
which the flux differences are expressed In terms 
of the mean value of the Jacobian matrix. The 
Jacobian matrix (1x1 In scalar case) Is required 
to satisfy 


(D f(u) - f(v) = A(u,v) (u-v) 
(11) A(u,u) = f u (u) = A(u) 


(2.14a) 

(2.14b) 


(111) A(u,v) has real eigenvalues and a (2.14c) 
complete set of eigenvectors 

Roe 3 has constructed a linearization of the form 
Eq. (2.14a) having these properties. Harten and 
Lax 7 show that such a Roe-type linearization 
exists as long as there Is an entropy function. 
Hence, from Eqs. (2.14) we can write 


ZffJ . f + ( Uj ) 


" a+( Y Yi )(u j - Yi> 


= a (u,,u, , ) A u 


J 


(2.15a) 


A + f' - r< Vl ) 


f '(Uj) 


j. U j +1 )(Uj +1 - Uj) 


< a Vx - 


<pt a 


j “j+1/2* f J 


+ Yl a j+1 /2 A f j 


AX 


(2.17) 


Hence we find 


,n+l 


= u" - x(A"f* + a + fT) n 


'](' - Xa j^1/2) a+f j 


- 0.5 xa 

-VI (• 

After rearranging, Eq. (2.18) becomes 


(1 + Xa j+1/2) A f j 


(2.18) 


Uj +1 = Uj - X 1 1 + 0.5 


+ + 
A f. 


♦ x h 



(1 ‘ Xa j +1 /2 )<p j 


A U ' 


-,n 


- (1 + xa 


A~f~ 



(2.19) 


a~(u 

= a"(u y u j+1 )A T Uj (2.15b) 

where a + (uj,uj_-j) >0, a-(uj,Uj + -|) <0. For 

scalar equation Eq. (2.1), the property, Eq. (2.14c) 
Is clearly satisfied. Thus the TVD conditions 
Eq. (2.12) can be ensured by this construction of 
a + and a-. This concludes the TVD consider- 
ation of the first-order scheme Eq. (2.8). 

Since first-order schemes contain too much 
dissipation, often obscuring the physics, we wish 
to convert them to TVD-satlsfylng, second-order 
accurate (both In space and time) schemes by add- 
ing limited antldlf f us 1 ve terms. Let us consider 
the one-step, Lax-Wendroff scheme 

Uj n+1 = uj - atfj * 0.5 at 2 (af x )JJ ♦ o(at 3 ) 

( 2 . 16 ) 

The third term on the RHS of Eq. (2.16) Is of 
second-order accuracy and Is considered as a part 
of the antldlf fuslve terms. We now can apply 
limiters to yield: 


The underlying Idea of the above arrangement Is to 
separately group " + " and fluxes. This turns 
out to be an Important step for the TVD proof. 

Now let 


L J-1 



j-1 /2' 


*4 - <P* 


(2.20a) 



1 - 0.5(1 + Xa^ +l /2^ 


1 + Xa 1-l/2 A f j - 



*>1 1 + xa 


J+l/2 A f j 


(2.20b) 


This provides us with an obvious choice for the 
functional relation for the limiter function, l.e., 


- ^ r J> 


(2.21a) 


where 


+ a 6^(1 - Xa j_l/2^ A Y 1 + Xa .1+1/2 ) 

r* = . . , r. = 


j a f j(l - xa J+1 /2 ) J a f j (1 + xa j-l/2 ) 

(2.21b) 


where 


Hence Eq. (2.20) becomes 


Combining Eqs. (2.23c) to (2.25) leads to 


+ x j 

I / + \ 

(* f j\\ 

1 + - xa j-i/2)(7 - V l]\ 

V VI 

1 'Vi /) 


(2.22a) 


0 < 



(2.26) 


Next, we Illustrate the entropy condition 
with a specific example. Let f(u) = u 2 /2, hence 
a = u. An obvious splitting of a gives 



).5^1 


+ \a 


3 + 1/27 




(2.22b) 


We note that Eqs. (2.19) to (2.21) are Identical 
In form to that obtained by Sweby*. Differences 
only appear In the detail; here we use the whole a, 
rather than a+ and a* In (r+, C + ) and (r“, C“). 

Since A~f + / A~u > 0 and (-A + f~/A + u) > 0, 

Eqs. (2.12a) and (2.12b) are satisfied by 
requl ring 


+ 



and Eq. (2.12c) leads to the condition 


(2.23a) 


(u ± |u|)/2 


(2.27) 


and 


f 1 = <T\i/2 (2.28) 

Note that since f Is a homogeneous function of 
degree 2, It Is easy to see that f 1 have con- 
tinuous first derivatives. Let the Initial 
condition 


u. x < o 

L t 


u o 


R. 


x > o 


(2.29a) 


and 


f(U L ) = f(U R ) 


(2.29b) 


The exact solutions are shown In Figs. 1. 



+ . _ + _ — . 

— X < 1 (2.23b) 

& u 

or 


\ | a | < 1 (2.23c) 

Note that in Eq. (2.23c), the TVD condition 
Eq. (2.12c) puts the same limit on the CFL number 
as the stability condition does. We note that 
Sweby 5 found: (1) the CFL number must be 

reduced to achieve the TVD conditions and (2) 

$ < 2. The differences between the present 
results and Sweby's are obviously due to the dif- 
ferences In formulation. 

We turn now to the determination of the 
bounds for the limiter functions. Allowing max- 
imum dlffusivlty to occur only In the first-order 
upwlndlng scheme requires 

q) 1 > 0 for any r± (2.24) 


For all TVD schemes using limiter-type formula- 
tions, an additional constraint Is Imposed, 

<p + - = 0, r± < 0 (2.25) 


That Is, a first-order scheme reverts as flux dif- 
ferences at neighboring points are of opposite 
signs, l.e., at point of extrema In flux. Hence 
second-order accuracy Is lost locally. This 
appears a price that has to be paid to make the 
limiter-type scheme TVD. Recent development of 
the so-called EN0( Essential ly NonOscI 1 latory ) 
scheme 8 remedies this shortcoming, but allows 
some "controlled" oscillations. 


x x 

centered- expansion stationary shock 

ul = - Ur < 0 u l = - ur > 0 

Figure 1 Exact (similarity) solution of Burgers' 
equation with Initial conditions, Eqs. (2.29). 

In the case, ul < 0, since the expansion jump Is 
not physically admlssable, hence centered- expansion 
fans give rise to unsteady solutions. On the other 
hand, the case ul > 0 allows a stationary shock 
solution for t > 0. It Is noted that Roe's orig- 
inal f lux-difference- spl It scheme and Murman- Cole ‘ s 
scheme for transonic potential equation admit 
expansion shocks. (See e.g., Harten et al. 9 ) 

In what follows, we see the performance of 
the first-order scheme Eq. (2.8) along with the 
splitting Eqs. (2.27) and (2.28). For calculation 
purpose, we assign 

( j < -1 

u . ( t - 0) = j ’ (2.30) 

J u R , j > o 

and take time step, X|u L | = 1.0. 

(1) Expansion, U|_ = - ur < 0 
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(1) Expansion, 

U L 

= - U R < 0 



time step u 

u 

-2 U -1 

u 0 

u 

At U L 

u 

L V 2 

Ur/ 2 

U R U | 

2At U L 

5u 

l /8 3 u l /8 

3Ur/8 

5 V 8 u . 

3flt 89 u l /126 

) u 

l /2 39 u l /128 39 u r /1 28 

Ur/2 89u ( 

Clearly the scheme 

does not admit an expansion 

shock. The expansions are centered about the 

midpoint between 

1 Initial jumps, x = -0. 

5ax, and 

spread over 2n 

grid points 

at t = nAt. 

(11) Shock, ul 

= - 

UR > 0 



time step 

U -2 

! U -1 

U o 

u i 

At 

U L 

0.5u l 

0.5u r 

U R 

2At 

U L 

0.75u l 

0.75u r 

U R 

3At 

U L 

0.688u l 

0 . 688u r 

U R 

4At 

U L 

0.71 5u l 

0 .71 5u r 

U R 

5At 

U L 

0.703u l 

0 . 703u r 

U R 

6At 

• 

• 

U L 

0.708u l 

0.708u r 

U R 

• 

OO 

U L 

0.707u l 

0.707ur 

UR 


We see that the shock wave Is stationary and Is 
represented with two Interior points. In general, 
both Gudonov's and Roe's schemes allow a station- 
ary shock with one Interior state while the 
Engqulst-Osher scheme admits two states. (See 
also van Leer 10 ) Since Roe's scheme also admits 
expansion shocks, the Gudonov scheme stands as the 
best for the representation of a stationary shock 
connecting the states Eq. (2.30), at least for 
f(u) = u2/2. 


For system of equations, the split-flux 
method Is shown to satisfy the entropy Inequality 9 , 
and Is thereby capable of selecting physically 
admlssable solutions. A contact discontinuity, 
however, will be smeared over a large extent, as 
seen later In Section 4. Harten 2 suggests use 
of artificial compression to further sharpen the 
contact discontinuity. However some care must be 
taken In a region of expansion so that entropy 
condition Is not violated, see Yee et al. 11 


We now express the present TVO scheme In 
terms of numerical fluxes 7j+-|/2 In the equation 

U j 1 = U j ‘ X( Vl/2 ' Vl/2* (2 ' 31) 


Rewriting Eq. (2.8a), the first-order upwind 
scheme yields 


J+l/2 




Vi + f j - 


2 , f 


f + - f" 


(2.32) 


We note that Eq. (2.32) Is not the Courant- 
Isacsson-Rees scheme 


’ j +1 /2 


■(' 


J+i 


f J - 


>1/2 




(2.33) 


where | a j+ 1/2 | 1s a mean value given by 

A + fj/A + Uj. Hence the case described by Eq. (2.29) 

will give rise to a^ ±]/2 = 0 and result In a 

stationary expansion shock, violating the entropy 
condition. 

After some algebraic manipulation, the 
second-order scheme Eq. (2.19) gives 




.1/2 * ( f j+l + f j - 4+f j * °M/2 4 * f j + >1/2 

(2.34) 


where 


J+l/2 = j^ 1 ’ Xa j+l/2^ " v j+l (• * Xa jfl/2^ 2 

( 2 . 

^ - Xa j^l/2) + ^l( ] + Xa J.l/2^ 


j +1 /2 


(2.35a) 
2 

(2.35b) 


The last two terms In Eq. (2.34) obviously contri- 
bute to second-order accuracy and are the antidif- 
fusion terms that also maintain TVO conditions. 


Remark 1. - Note that 1 + ^ a j+l/2 > 0 for 
stability, hence Gj + -|/ 2 satisfies the 
Inequality 


1 > °j + i/2 > °* 1f 0 < < 1 (2.36) 

and Is of opposite sign to the third term on RHS 
of Eq. (2.34). Thus, It is clear that the present 
and other TVD schemes share the common feature In 
which the higher-order TVO procedure amounts to 
adding more sophisticated antidiffusion terms to 
a basic lower-order scheme which Itself Is a TVO 
scheme. 


Remark 2. - On the other hand, It Is more 
advantageous to think of It as adding diffusive 
terms to an underlying higher-order, non-TVO 
scheme so that TVD conditions are enforced. The 
one-step Lax-Wendroff scheme defines 

? JO/2'( f M ")- »).,//' j)/ 2 < 2 ”> 

In Eq. (2.31). Rewriting Eq. (2.34) gives 

>1/2 ■ f J+l‘/2 * K 1 - - ( a j.l/2 * Xa j.l/ 2 ) 4 ' f j |/ 2 

(2.38) 

where It can be shown 


for 0 < <p~ < 1 . The last bracket In Eq. (2.38) 
as substituted In fj + i/ 2 - fj-1/2 seen to yield 
a nonlinear second-order central difference. 
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(3.6a) 


Remark 3. - Note that <Vl/2 A+f j and 
(1 - • 3 j +1 / 2 ) A+f j 1n Ec l' ( 2 - 34 ) correspond to 
(g^ + g^ +1 ) and Qj + -| /2 a+u j 1n Harten’s second- 
order scheme. 2 - 11 


A = A* + A , a^ = a^ + a^ 

A 4 = (A ± |A|)/2, jA^jj = ja, j 
Consequently A Is split as 


(3.6b) 


Remark 4. - For smooth solutions of Eq. (2.1), 

we have r 4 = 1 + 0(|/Tu|) and <p 4 = 1 + 0(|A ± u|). 
It Is easy to show that 


A = A + + A 


(3.7a) 


A 4 = SaV 1 2 


(3.7b) 


♦ 0(I>1 “ |:| > (2 - 39) 

hence Eqs. (2.19) or (2.38) Is second- order 
accurate. 


S yste m of Conserva tion Laws- 
O ne- Dimensio nal Equations 

In this section we describe how to extend the 
procedure developed for the scalar equation, 
namely upwlndlng combined with the notion of lim- 
iters, to systems of conservation laws. To show 
the extension, we will exploit the fact that In 
one space dimension, the system Is hyperbolic and 
hence Is diagonal 1 zable. As the system Is 
decoupled, the mathematical procedure essentially 
follows what has been described In Section 2. But 
It Is less clear as to how to define the argument 
r for the limiter tp, which Is a key element In 
the limiter-type schemes. 

Let us consider the one-dimensional system 


U t + F X (U) =0, t > 0 
U(x,0) = U Q (x) 


(3.1) 


Here U and F(U) are column vectors of m com- 
ponents, the flux vector has a Jacobian matrix 


We note that to arrive at the split in matrix A, 
Eq. (3.7) the only condition used is that A is 
diagonal 1 zable. We have as yet made no connection 
with the flux vector F. 

Now let A be constant, hence S and A 
are constant. From Eq. (3.2) we have 

F - AU (3.8) 

for any U. With the aid of Eq. (3.7), Eq. (3.8) 
gives a split In F 


F = F + + F~ (3.9a) 

F 4 = A 4 U (3.9b) 

Here we remark that as a result of assumption (2), 
F and F- are homogeneous functions of degree 
one. On the other hand, for the nonlinear system, 
specifically the Euler equations, Steger and 
Warming^ 2 need to rely upon the homogeneous 
property of F to accomplish the flux splitting. 

Employing the assumption (1), we define a 
vector W by 

$ _1 U = W, W = (w l : l - 1 ,m) (3.10) 

Hence 


A(U) 


3FLU1 

au 


Eq. (3.1 ) is rewritten as 

U t + A(U)U x * 0 


(3.2) 


(3.3) 


To complete the extension and the proof of TVD 
for the system, we need to make the following 
assumptions : 2 

(1) A has real eigenvalues and a complete set 
of eigenvectors, hence A is diagonal 1 zable , 

(2) A is constant. 


S ] F = AW 

(3.11a) 

-1 + + 

» F“ = A~W 

(3.11b) 


Following the same procedure described In 
Section 2, the spatial derivatives which appeared 
in the Lax-Wendroff scheme are now approximated by 

. (** f i • .ill) . .feL‘l F L-.ridl*!iI 

X = AX ° 2AX 

(3.12) 


and 


Let S be a matrix, the columns of which are the 
right elgn evectors of A. Then 

S _1 AS = A (3.4) 

where A Is a diagonal matrix 

= a i 4 ij> 1 £ 1 J £ m (3.5) 

Since A Is real, we can write 


< AF x>x - * 


_ A - / W( <t> .1 A * F .i * »3 +/ F j) 


AX 


(3.13) 


Here cp- are scalar functions, to be determined 
later. Then we have 


6 



Cj = (-xa ) 


,,n+l , ,n . - r + ,V\n 

u j = u j ■ x(a F j + a F j ) 

- 0.5 XA _ j*J(I - XA J+1/2 )A + Fj 

- ‘ | " < 3,, > 

Letting cp 1 = 0 reduces Eq. (3.14) to the Euler 
explicit, first-order, upwind scheme 

Uj +1 = Uj - \(a"Fj + A + Fj ) (3.15) 

In the following proof of TVD conditions for 
a system, Eq. (3.14) is considered and we will 
treat Eq. (3.15) as a subset of Eq. (3.14), 
although the second-order accurate TVD scheme is 
really built on the first-order scheme, a point 
brought out in the previous section. Premultiply- 
ing Eq. (3.14) with $ - and using Eqs. (3.10) 
to (3.11), one gets 


w n+1 w n 

W J ' W 3 


- 0.5 XA 


- \(a + a'Wj + a aVj) n 

|tpj(I - XA) A + A + W 

(3.15) 


j - <pj +1 (I - XA)A aVj " 


Since A- are diagonal matrices and <p j * s are 
scalar functions, Eq. (3.16) is a m-component 
system of decoupled equations. 


(WjV +1 = (w}) n - X(a; 4 -«; * a-A + WjV 
- 0.5 XA |(pj(l - xa^Ja^A^Wj 

‘ x a t )a i A+w j | n » 1 1 a < m (3.17) 


This is identical to the constant-coefficient 
version of the scheme of Eq._(2.19) for the non- 
linear scalar equation if A+f 1 = a ± A + u. Proof of 
TVD however departs here from that in Section 2. 
The most important thing to keep in mind is: the 

transformation (decoupling) is employed merely as 
a mechanism to ensure TVD, the ultimate goal is 
treating the system. Therefore, in the process of 
proving TVD conditions for Eq. (3.17), definitions 
regarding r and (p(r) must not hinder the 
transformation back to the coupled system of 
Eq. (3.14). 


Dropping superscript "fc," Eq. (3.17) is 
rewritted as 


with 


n+1 n _+ - n r - + n , os 

w j = - c j-i a w j + c j a w j (3 - 18) 


c j_l 55 Xa+ h + 0*5(1 - Xa) 


+ + 

a A W ,1 ♦ ♦ 

~ *1 " *J-1 

a A 


(3.19a) 


1 - 0.5(1 + xa) 


a A w. 

" 7aV, "3 


(3.19b) 


Here the trick is: associate all term having "+" 

eigenvalues with <p + and those of H - H eigen- 
values with q>~ . Also, recognizing the fact 
that the transformation is possible only for those 
w connected with a^, it becomes clear that we 
must redefine r for the system. An obvious 
choice is 


r . L£’i 
) ■ •-» 


a A w 


3 


- +■ 

a A w. 


(3.20a) 


(3.20b) 


a A w 


3 


Eqs. (3.19) are rewritten as 


cj_, = (xa + ) 


Cj = (-xa ) 


1 + 0.5(1 - xa) 


V r j 


- <p 


1 - 0.5(1 *■ xa) (<P 


'J+l 




r J> 


(3.21a) 


(3.21b) 


Harten's TVD conditions Eqs. (2.12) are satisfied if 

= $ < 1 (3.22) 


+ 

!l ± 

r ± ' 
r i 


and 


x | a | <1 


(3.23) 


as Eq. (2.23). This completes the TVD proof for 
each component of the decoupled system. 

If a solution procedure for Eq. (3.1) is con- 
structed for the decoupled system Eq. (3.17), we 
see that the definitions of Eq. (3.20) require 
knowledge of the eigenvalues and eigenvectors of 
A (see Eqs. (3.4), (3.5) and (3.10)). This causes 
concerns not only in numerical operation counts 
but also in extention to multidimensions, because 
decoupling systems in multidimensions may be an 
extremely complex task if possible at all. Hence, 
the goal is to express r± in terms of simple 
variables e.g., U or primitive variables, or pre- 
ferably fluxes. The latter has the attractive 
property being able to naturally reduce to the 
scalar case. 

We turn now to show a rational procedure to 
arrive at definitions of r 1 for systems of 
equations. Putting the superscript "fi." back in 
Eq. (3.20a), (r- Is done similarly) but dropping 
subscript "j", we have 
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a*A w l = r*a*A + w l , 1 < 8. < m 


or In matrix form 


aVw = rVa + W 


(3.24a) 


(3.24b) 


where R + Is a diagonal matrix whose diagonal 
elements are { r { : SI □ 1 , m} . Premultiplying 

both sides with S brings us back to the coupled 
system, 


aYu = rVa + u 


(3.25) 


With the aid of Eq. (3.9), we can relate fluxes 
to r + 

(3.26a) 


A'F + = rVf + 


or 


A F Sl = r i A+F I * 1 1 i < m 


(3.26b) 


Since the signs of the flux differences are key 
factor In determining the limiter function q>, use 
of norms must be ruled out. However, some kind of 
sum Is needed In order to Include Influences from 
each component. Several possibilities Immediately 
follow from Eqs. (3.26) . 


(a) I>+A + F + -Za’fJ 

SL * * S i * 


(b)2>'(aY) 2 


(O r[ 

l * 11 


1 + 
4 F * 


a + f ; = E a'fj 

si 1 


A ' F I 


(3.27a) 

(3.27b) 

(3.27c) 


The last choice mimics the L 2 -norm but retains 
the sign. We see that Eqs. (3.27) offer the pos- 
sibility of using different weights for different 
components of F. Let all components be of equal 
weight, then from Eqs. (3.27) we can define: 


(a) r + = = < a~ F + . 1 > j <aV, 1 

(3 

* , < 4 V. 4 V >/< A + F + , »V 


(3.28a) 


(b) r 


(3.28b) 


(c) r + = Z| A " F i| - < | A ‘ F+ | • A~ F + > / < | 4 * F *|’ A+F+ > 

(3.28c) 

(d) let 


+ 

r l - 

r + of 

Eq. (3.28c) 


and 





[°- 

rj < 0 


+• ' 
r = j 

) 

1 . f 

1/2 + 

( 3.28d) 

I 

' (rp 

, r > 0 



where <,> denotes Inner product among all m com- 
ponents. Clearly, the first choice Is simplest. 

Now, since <p- are scalar functions, what Is 
used for the scalar equation carries over automat- 
ically to systems. This completes the construc- 
tion of TVD scheme Eq. (3.14) for systems of con- 
servation laws. 

Remark 1. - As In the scalar case (Section 2), 
we denote In Eq. (3.14) 


+ + 

= <P(I%) 




(3.29) 


Remark 2. - In the course of developing the 
present scheme, the assumption of A being con- 
stant Is necessary. For the nonlinear case, a 
formal generalization Is made. Thus, the same 
Integration scheme of Eq. (3.14) still applies 
since It Is the flux differences that are used 
and no transformation Is needed for evaluating r 
and <p. Hence, the present scheme Is not limited 
to flux-vector splittings (for example, Steger and 
Warming 12 and van Leer 13 ), the Roe- type flux 
difference splitting 3 Is equally applicable. 

Remark 3. - We avoid matrix-vector multipli- 
cations as needed In Harten‘s scheme or a Rlemann- 
type solver (See e.g., Osher and Chakravarthy 4 , 
Roe 3 ), but require the calculation of " + " and 
fluxes. The construction of "+" and flux dif- 
ferences Is as easy as that of A-F , the differ- 
ences In whole flux. Also, the calculation 
required of r 1 and cp± Is minimal once A ± F i 
are formed; the amount of operations In Eq. (3.14) 
should be nearly twice as much as the 
convenntlonal 

non-lVD, l.ax-wendrof f scheme(w1thout adding addi- 
tional dissipation). Also the Increase In storage 
varies with schemes, and Is minimal for explicit 
schemes . 

^ Remark 4. - We can write the numerical fluxes 
Fj±l/2 f° r the system In the same way as the 
scalar case (Section 2). If fact, one merely 
replaces In Eqs. (2.32) to (2.35) the scalar flux 
with flux-vectors and a with the Jacobian matrix 
A. 

Gener a_l Iza 1 1 o n_ t o_ Eul er _ Eg uat Ion s - 
Multidimensions jand^ Steady Equations 

In this section we apply what has been 
described In Section 3 toward the solutions of 
Euler equations In one- and two-space dimensions. 
To describe the extension, we consider the 
one-dimensional equation 


U t + F X (U) = 0 


where 


U = 


p 


pU 

pU 

, F = 

pu 2 + p 

pE 


(pE + p)u 


(Y - 1) P(E - 0.5 U 2) 


(4.1a) 


(4.1b) 


(4.1c) 


To accomplish the differencing of F + and 
F" In Eq. (3.14), we choose to use flux-vector 
splittings given by Steger-Warmlng 1 2 and van Leer 13 . 
They are written, for completeness, as follows: 
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(a) Steger- Warming 

Let the eigenvalues of A be 

2 

(a r a 2 ,a 3 ) = (u,u + c, u - c), c = yp/ P (4.2) 


and 


+ 



(*« * 


1 < JL < 3 


Since the flux vector In Eq. (4.1) Is a homogeneous 
function of degree one, we have 

F = AU, A = 8F/3U (4.3) 


Hence, the flux Is split as 


F 1 = A 4 U = §- 
2y 


2(y - 1)3^ + a 2 + a" 

. . + + + 

2(y - l)a 1 'a 1 + a~a 2 + a“a 3 

4 ? / 4 ? 4 ?\ 

( y - l)a"a-j + 0.5 la“a 2 + a"a 3 1 + w 


(4.4) 


where w = (3 - Y)(a 2 + a*)c/2(Y - 1)* 


Figures 2 show the comparisons of calculated and 
exact solutions of pressure, density, Internal 
energy, and velocity; the Initial condition Is 
also Included. The Steger-Warmlng splitting was 
used; the CFL number was 0.95 and the domain 
divided equally Into 200 Intervals. The limiter 
function was 


( m1n(l ,r) , r > 0 
(0, r < 


where r was defined by Eq. (3.28c). It Is evi- 
dent that the contact discontinuity was smeared 
due to excessive dlffuslvlty. Harten 2 reported 
substantial Improvement In the resolution of the 
contact discontinuity through the addition of arti- 
ficial compression. Here we show the results of 
using different definitions given In Eqs. (3.28). 
(Note: the definition for r-, although not given 

here, can be easily gotten from Eq. (3.20b) and by 
following steps Eqs. (3.21) to (3.28)). The defi- 
nition of Eq. (3.28a), the results of which are 
shown in Figs. 3, clearly gave the best results 
while the results corresponding to Eq. (3.28b) 
were the poorest, displaying as evident In Fig. 4 
slight oscillations near the contact discontinuity. 
However, the shock wave and expansion waves were 
predicted very well for the nonlinear system In 
all cases. Next Roe's "superbee" 14 was used. 


(b) van Leer 

To circumvent the discontinuity of first der- 
ivatives of F 1 at transition states, u = 0, 
or u = ±c, van Leer constructs a flux-vector 
splitting consisting of a polynomial of degree 
two, thereby ensuring continuous first-order 
derivatives . 


For supersonic flow, |u| > c 
F + = F , u > 0 

F' = F , u < 0 

For subsonic flow, |u| < c, let 
F^ = ± p(u i c) 2 /4c, 


(4.5a) 


4 



= F* [( Y - l)u 1 2 c]/y 
f] [(Y - l)u i 2 c] 2 /2(y 2 


(4.5b) 


D 


Both splittings were used in the numerical 
experiments to determine their appl 1 cahl 1 1 ty to 
the Dresent TVD scheme. In the following, we dis- 
cuss the specific problems studied. 


I max jm1n(2r,l), m1n(r,2)j, r > 0 

(4.8) 

0, r < 0 

It Is seen In Fig. 5 that the contact discontinuity 
was further improved. But slight oscillations 
occurred in the expansion region, possibly due to 
the fact of using $ < 2 which violates the TVD 
conditions In the present formulation. We also 
devised a mechanism by which a dramatic improve- 
ment was observed at the contact discontinuity; 
more analysis is underway and will be reported 
elsewhere. 

Next we tested a case Involving a much 
stronger shock wave with the Initial data: 

((500, 400, 0.0), x < o 
(P,p,u) = < (4.9) 

((1.0, 1.0, 0.0), x > o 

Figures 6 show the results calculated using 
Fq. (3.28a); the shock wave as well as the expan- 
sion waves were accurately predicted and monoto- 
nicity was satisfied. However, the extent of 
smearing at the contact discontinuity was not 
acceptable. Hence, calculations with finer meshes 
of 500 equal spaclngs were made, as seen In 
Figs. 7 much better results were obtained. 

Qua s 1_ One -_D 1 me n s j on a 1_ Hozzje froblem 


0 ne- D i men s Ion aj .Shoe k Tube f rob Tern 


The flows are described by 


lo test the performance of the second order 
accurate scheme Eq. (3.14), first we consider the 
Sod problem with the initial condition. 


(P^,u) 


( 10 . 0 , 8 . 0 , 0 . 0 ), 

( 1 . 0 , 1 . 0 , 0 . 0 ), 


x < o 
X > 0 


(4.6) 


U t + F x (U) = H(U), H( U) 


o 

p dS/dx 
o J 


(4.10) 


where S Is the cross-sectional area of the 
nozzle. Two nozzle shapes were considered: 
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S^x) * 1 .398 + 0.347 tanh (0.8x - 4), 0 < x < 10 

(4.11a) 


and 


( 1 . 75 - 0.75 COS ( 0 . 2x - l)ir, 0 < x < 5 
S ? (x) = ] 

* (1.25 - 0.25 cos (0.2x -l)ir, 5 < x < 10 

(4.11b) 

Uniform grids of 100 Intervals were used In 
both nozzles. In a divergent channel where the 
Inflow Mach number was 1.26 and the exit pressure 
ratio was 0.746, the van Leer splitting Eq. (4.5) 
edged Steger-Warmlng ' s In resolving the shock 
wave. Fig. 8. This Is related to the continuous 
transition of fluxes splitting through the sonic 
state inside the "numerical" shock wave. It Is 
more evident In the convergent-divergent channel 
where the Inflow Mach number and the exit pressure 
ratio are 0.2395 and 0.84 respectively. Figure 9 
displays clearly the "kink" at the sonic throat In 
Sterger-Warmlng ' s splitting. The van Leer split- 
ting performed admirably well at the throat 
although the exact solution has a discontinuity in 
first derivative resulting from a jump In $ 2 ". 

We also applied the same Idea for construct- 
ing the above-described TVD scheme, to the 
MacCormack explicit method. The purpose was to 
see whether and how well the present Idea worked 
In a scheme where each of the predictor- corrector 
steps utilizes one-sided (but not upwind) differ- 
ences. The calculated results, although not 
shown, were essentially Identical In both nozzle 
flows . 

Tw o- D1 mensionaj JP_roJ)lern 

A completely different approach toward the 
solution was taken: (1) Steady- state equations 

were used, 1 .e. 


F X (U) + G y ( U ) = 0 (4.12) 


Here the type of the equations may be no longer 
hyperbolic and hence the notion of TVD Is unclear; 
(2) We used Newton's method to linearize the sys- 
tem Eq. (4.12) and obtained an Implicit Integra- 
tion procedure, as opposed to the Lax-Wendroff 
type; And finally, (3) The basic underlying scheme 
used second-order upwind differencing and was con- 
verted to a form similar to Eq. (3.12). For 
example, the derivative of F Is approximated by 



+ A F, 


+ 


1 

2 



- ,-A“F 



AX 


(4.13) 


Without mathematical guidance, we simply 
carry over the preceding formulation. However, It 
Is assumed that since the time-asymptote steady 
solution to the unsteady equations can preserve 
monotonicity, the solution to the steady equations 
must also have the Identical property (uniqueness 
of solutions is assumed). 


In one-dimensional problems, the shock wave 
Is aligned with a grid line and can be represented 


with a sharp profile by TVD schemes. In this sec- 
tion we wish to test our scheme on two-dimensional 
problems. In particular where the shock waves cut 
across grid lines. We consider a regular reflec- 
tion of an oblique shock wave from a solid surface. 
Two Inflow Mach numbers were tested to see the 
effect of shock strength on the performance. The 
Inflow conditions were fully specified with free- 
stream values and the conditions at the top bound- 
ary were set to satisfy the shock-jump relations 
with a specified shock angle. The variables at 
the outflow boundary were extrapolated linearly. 

At the solid wall, we let v and the gradient of 
the other variables vanish. (See also ref. 11). 

The first case had an Incoming Mach number of 
2.9 and a shock angle of 29.0 degrees; the computa- 
tional domain contained 60x20 meshes equally divided 
in a domain 0.0 < X < 4.0, and 0.0 < Y < 1.0. 

Figures 10 show the comparisons of pressure from 
different calculations with exact solutions at 

Y = 0 and 0.5, together with pressure contours. 

The second-order upwind scheme with no added 
artificial-damping displayed over-expansion and 
over-compression just upstream and downstream of 
the shock, thereby causing slight Irregularities 
In contours near the shock wave. All calculations 
with the present procedure did predict monotonic 
transition through shock waves, but In most cases 
the shock waves were smeared to a much larger 
extent than In the one-dimensional problems. The 
definition of Eq. (3.28a) gave slightly better 
representation of shock waves than Eq. (3.28c). 

The best came with using van Leer's splitting and 
Eq. (3.28a), showing much Improved shock-wave 
resolutions, perhaps as good as from the one- 
dimensional cases. This might look surprising 
since in the supersonic flows both splittings 
should give the Identical results. However, this 
Is true only for the x-flux F, but the y-flux G 
Is split differently as the transverse velocity 
component Is subsonic. 

Next let us consider a strong shock wave 
resulting In an overall pressure rise of a factor 
of 153; the Incoming Mach number was 10.0 and the 
shock angle remained the same as above. We divided 
the domain 0.0 < X < 5.5, and 0.0 < Y < 1.0 Into 
82x20 meshes with same AX as above. The van 
Leer splitting was used along with definition 
Eq. (3.28a). The pressure distributions and con- 
tours are given In Figs 11. The pressure at 

Y = 0.5 was remarkably well-behaved. However, 
the monotonicity seemed lost In the surface pres- 
sure, possibly due to extrapolation from the 
interior solutions. Roe's "superbee" displayed 
slightly sharper representation of shock waves 
than the case $ < 1.0, but also developed a 
stronger overshoot In the surface pressure just 
behind the shock wave. 

C oncluding Remarks 

In the present paper we showed a simple and 
general procedure for constructing a TVD differ- 
encing scheme In which second-order accuracy was 
maintained, except in regions having extrema of 
fluxes. We showed that any non TVD dlf f erenci ngs 
could be converted by writing It as a first-order 
upwind scheme plus remaining higher-order fluxes 
differences upon which limiting functions were 
applied so as to satisfy Harten’s TVD conditions 
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In the nonlinear scalar equation and constant- 
coefficient hyperbolic system of equations. A 
mathematical proof was given and definitions of 
arguments for the limiter functions derived for 
the system. Generalizations to the nonlinear 
Euler equations were carried out. Numerical tests 
were conducted to see the performance of the pro- 
posed scheme; questions regarding limiters were 
Investigated. We found the definition of the 
argument appearing In a limiter was quite cri- 
tical In getting good resolution for discontinu- 
ities, but otherwise all definitions seemed to 
yield monotone behavior. Further Improvement In 
representing a contact discontinuity was also 
found possible. It has become clear to us that 
more analyses can be done In the area of limiters 
to achieve substantial Improvements. We shall 
report this In a forthcoming paper 15 . It was 
also demonstrated that for steady-state calcula- 
tions the van Leer splitting was In all cases more 
accurate than the Steger-Warmlng splitting. 
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FIGURE 2. - ONE-DIMENSIONAL SHOCK TUBE PROBLEM. INITIAL CONDITION AS GIVEN 
IN EQ. (4.6) . RESULTS FROM THE SOLUTION USING EQS. (3.28c) AND (4.7). 
AX = 0.05, CFL = 0.95. 
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FIGURE 3. - ONE-DIMENSIONAL SHOCK TUBE PROBLEM, INITIAL CONDITION AS GIVEN 
IN EQ. (4.6). RESULTS FROM THE SOLUTION USING EOS. (3.28a) AND (4.7), 

AX = 0.05, CFL = 0.95. 


PRESSURE DENSITY 



PRESSURE DENSITY 




PRESSURE DENSITY 

>0375 500 0 80 160 2*0 320 *00 



PRESSURE DENSITY 





FIGURE 7. - ONE-DIMENSIONAL SHOCK TUBE PROBLEM, INITIAL CONDITION AS GIVEN 
IN EQ. (i). 9). RESULTS FROM THE SOLUTION USING EQS. (328a) AND (9.7), 

AX = 0.02, CFL = O.A75. 
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(A) STEGER-WARMING SPLITTING, (B) VAN LEER SPLITTING, 

EQS. (4.4). EQS. (4.5). 


FIGURE 8. - STATIC PRESSURE DISTRIBUTION IN THE DIVERGENT NOZZLE, EQ. 
(4.11a), Moo = 1.26, RATIO OF STATIC TO TOTAL PRESSURE = 0.746. 



(A) STEGER-WARMING SPLITTING, (B) VAN LEER SPLITTING, 

EQS. (4.4). EQS. (4.5). 


FIGURE 9. - STATIC PRESSURE DISTRIBUTION IN THE CONVERGENT-DIVERGENT 
NOZZLE, EQ. (4.11b), M* = 0.2395, RATIO OF STATIC TO TOTAL PRES- 
SURE = 0.84. 



(A) SECOND-ORDER UPWIND SCHEME, SWS. 

FIGURE 10. - SHOCK REFLECTION PROBLEM, STATIC PRESSURE 
DISTRIBUTION AT y = 0 AND y = 0.5 AND PRESSURE CON- 
TOURS (MIN = 0.5, MAX = 4.0, INCREMENT = 0.007), 

Moo ~ 2.9, SHOCK ANGLE = 29 DEGREES. 


(B) EQS. 


FIGUR 





O 

O 

S-H SPLIT 
OX * 0.067 
♦ ♦♦♦ EXACT 
o c o • COMP 



J — ± 


X/L 


{ 




meaBmammimeam 



© 

© 

o 



S-H SPLIT 
OX = 0.067 
♦ ♦♦♦ EXACT 
o o o o COMP 


J — i — 1 — i — 1 — { 

X/L 



.28c) AND (A. 7), SWS. 
10. - CONTINUED. 





(C) EQS. (3.28a) AND (4.7), SWS. 
FIGURE 10. - CONTINUED. 





(D) EQS. (3.28a) AND (4.8). SWS. 
FIGURE 10. - CONTINUED. 




(E) EQS. (3.28a) AND (4.8), VLS. 
FIGURE 10. - CONCLUDED. 
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(A) EOS. (3.28a) AND (4.7), VLS. 

FIGURE 11. - SHOCK REFLECTION PROBLEM, STATIC PRESSURE 
DISTRIBUTION AT y = 0 AND y = 0.5 AND PRESSURE CON- 
TOURS (MIN = 0.5, MAX = 150.0, INCREMENT = 2.25), 

Mo, = 10.0 AND SHOCK ANGLE = 29 DEGREES. 


(B) EQS. (3.28A) AND (4.8), VLS. 
FIGURE 11. - CONCLUDED. 
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